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Abstract 

A new numerical method for an inverse problem for an elliptic equation with un- 
known potential is proposed. In this problem the point source is running along a 
straight line and the source-dependent Dirichlet boundary condition is measured as 
the data for the inverse problem. A rigorous convergence analysis shows that this 
method converges globally, provided that the so-called tail function is approximated 
well. This approximation is verified in numerical experiments, so as the global con- 
vergence. Applications to medical imaging, imaging of targets on battlefields and to 
electrical impedance tomography are discussed. 



1 Introduction 

The phenomenon of multiple local minima and ravines of least squares residual functionals 
represents the major obstacle for reliable numerical solutions of Multidimensional Coefficient 
Inverse Problems (MCIPs) for Partial Differential Equations (PDEs). We believe that be- 
cause of the applied nature of the discipline of Inverse Problems, the issue of addressing the 
problem of local minima has vital importance for this discipline. Indeed, any gradient-like 
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optimization method for such a functional would likely converge to a local minimum, which 
is located far from the correct solution. Furthermore, a global minimum, even a well pro- 
nounced one, is not necesseraly located close to the true solution, because of the ill-posed 
nature of MCIPs. Because of this, the vast majority of current numerical methods for MCIPs 
are locally convergent ones, like, for example Newton-like method, see, e.g., [1],[2],[4],[5] and 
many issues of Inverse Problems. That is, convergence of such a method to the true solution 
is rigorously guaranteed only if the initial guess is located sufficiently close to that solution. 
However, in the majority of applications such as e.g., medical and military ones, the media of 
interest is highly heterogeneous, which means that a good first guess is unknown. The latter 
naturally raises the question about the reliability of locally convergent numerical methods 
for those applications, and this question is well known to many practitioners working on 
computations of real world MCIPs. 

Thus, we are interested in the issue of globally convergent numerical methods for MCIPs. 
We call a numerical method globally convergent if the following two conditions are in place: 
(1) a rigorous convergence analysis ensures that this method leads to a good approximation 
of the true solution regardless on the availability of a first good guess, and (2) numerical 
experiments confirm the said convergence property. 

In this paper we present an "almost" globally convergent method for an MCIP for the 
equation 

A x m - a (x) u = -5 (x - xo) , x = (x, z) G R 2 , (1.1) 
lim u(x,xo) = 0. (1.2) 

Here x is the source position, and this position is running along a line to generate the data 
for the inverse problem. We use the word "almost" , because we rigorously prove global con- 
vergence only assuming that we know a good approximation for the so-called "tail function" , 
i.e., we assume that we know a good approximation of the second term of the asymptotic 
behavior of the function In [w (x, x )] for |x | — > oo. Since this approximation is unknown 
analytically, we have decided to use a heuristic iterative "accelarator" for convergence of 
tails and to confirm the desired convergence numerically. Assuming that our tail function is 
close to the correct one, we prove a global convergence result, which does not rely on a good 
first guess for the solution. This is why we call our method "globally accelerated". From 
the numerical standpoint, another advantage for using the accelerator is that it gives us an 
approximation for the tail, which seems to be rather close to the actual tail and we observe 
this numerically. The only drawback is that we cannot establish this rigorously. 

We assume throughout this paper that the function a G C a (IR 2 ) , a > const. > and 
o(x) = k 2 = const. > for x G R 2 \f2, where a G (0, 1) . The classic theory implies that there 
exists unique solution of the problem (1.1), (1.2) such that u G C 2+a (|x — x | > e) , Ve > 0. 

The first generation of globally convergent numerical methods has started from the so- 
called convexification algorithm [6]. The convexification is using the projection technique 
with respect to all variables, except of one, and a stable layer stripping procedure with respect 
to the latter variable. While the work [6] is concerned with time/frequency dependent 
data, our publication [8] is applying the convexification to the case of the running source, 
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which is the same as one in this paper. In the mathematics literature some other numerical 
techniques providing the global convergence (see, e.g., [15-20])) are available. Their numerical 
implementations can be found in [21-24]. 

Recently the development of the second generation of globally convergent numerical 
methods was initiated in [3]. The idea of [3] was originated from our earlier publication [9]. 
The concept of [3] also overlaps in part with the scheme of the current paper. However, 
unlike our current case, the time dependent data resulted from a single measurement are 
considered in [3]. Laplace transform of either hyperbolic or parabolic equation of [3] leads 
to the equation Aw — s 2 c(x)u> — a(x)w = — 8 (x — x ) , where s > is the parameter of 
the Laplace transform and the source position is fixed. Compared with our case, the main 
advantage of this equation is that the asymptotic behavior of tails at s — > oo is known. 
Specifically, in the case when the coefficient c(x) is unknown, lim,,^^ (In w/s 2 ) = 0, and 
similarly when a(x) is unknown. Ultimately, the knowledge of these limits enables one to 
prove a global convergence theorem in [3] without a heuristic assumption. 

We are interested in the extension of the idea of [3] to the case of the running source 
instead of the changing time or frequency. In other words, we consider almost the same 
inverse problem as one in [8]. However, instead of the convexification of [8] we develop an 
analogue of the method [3]. A numerical method, similar with one of this publication, was 
published in our early work [9]. However, the treatment of tails in section 4 of [9] was 
different from one of our case, and that is why the global convergence property was not 
observed in [9]. We also refer to subsection 5.4 of [4] for another treatment of tails for a 
Newton-like locally convergent method for an MCIP with frequency dependent data. We 
now explain the underlying reason of our difficulties with the tail function from the physics 
standpoint. In the case of the time dependent data for a hyperbolic equation [3] the tail 
function is close to the so-called "first arrival wave" . It is well known that the first arrival 
signal is very informative one. However, it is unclear what the first arrival signal is in our 
case of the elliptic equation (1.1) with the running source. 

We now formulate our inverse problem. 

Inverse Problem. Denote x = (x, z) . Let Q C 1R 2 be a bounded domain and T = dQ. 
Let B be a constant. Suppose that in (1.1) x = (s,B) Q. Determine the coefficient a(x) 
for x G VL, assuming that the following function ip (x,s) is given 



where s is a sufficiently large number, s<sisa certain fixed number and {x = (s, B) , s > s}fl 



We consider the 2-D case for the sake of simplicity only for this complicated problem. 
Generalizations of our method on the 3-D case are feasible. The parameter count shows 
that the data </?(x,s) depends on two free parameters, so as the unknown coefficient a(x). 
Hence, this Inverse Problem is non-overdetermined. This inverse problem has applications in 
imaging using light propagation in a diffuse medium. This is the so-called continuous-wave 
(CW) light. In this case the coefficient a(x) is 




(1.3) 



dfl=0. 
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where /x' s (x) is the reduced scattering coefficient and /i (x) is the absorption coefficient of the 
medium [1] . The first example of this application is imaging of targets on battlefields covered 
by smog and flames using propagation of light originated by lasers. In this application 
the laser source should be moved along a line and the measurements of the output light 
should be performed at the boundary of the domain of interest. Interestingly, the diffuse- 
like propagation of light would be even helpful, because the direct light can miss the target. 
The second applied example is in imaging of human organs or small animals using near 
infrared light propagation. Note that this application is discussed in many publications, in 
which locally convergent numerical methods are developed, see, e.g., [1],[9],[13]. Also, the 
above Inverse Problem has applications in Electrical Impedance Tomography, in which case 
the original equation is V • (a (x) Vt> ) = — 5 (x — x ) and the standard change of variables 
u = V\fa reduces this equation to (1.1) assuming that a (x) = 1 in a neighborhood of the 
source position x . Here a (x) > const. > is the electric conductivity of the medium. 



2 Nonlinear Integral Differential Equation 

Since the function u is positive, the by the maximum principle we can consider the function 
v — lau. Since the source x = (s, B) ^ Vt, we obtain the following equation from (1.1) 

Av + \Vv\ 2 = a(x) in Q, (2.1) 

v (x, s) = ^ (x, s) , V (x, s) e r x (A, s) , (2.2) 

where tp 1 = Imp. To eliminate the unknown coefficient a(x) from equation (2.1), differentiate 
it with respect to s and let 

q (x, s) = d s v (x, s) . (2.3) 

Then 

s 

v (x, s) = - J q (x, r) dr + T(x), x eQ,s e (s, s] (2.4) 

s 

In (2.4) T(x) is the so-called "tail function". The exact expression for this function is 
of course T(x) = v (x, s) . We know only the first term of the asymptotic expansion of the 
function v (x, s) at s — > oo (below). As it was pointed out in Introduction, if we would know 
the second term also, as it is the case of the time dependent data of [3], then we would 
be better off approximating the tail function. However, the absence of the knowledge of 
this term significantly complicates the matter compared with [3]. Thus, we develop below 
a heuristic iterative procedure of an iterative approximation of the function T(x), with the 
aim of funding such an approximation T appr (x) that VT appr (x) w VT (x) . 

We obtain from (2.1)- (2.4) 

s 

Aq - 2Vg J Vqdr + 2VgVT = 0, (2.5) 
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g(x,s) = V(x,s),V(x,s) G T x (A, s) , 



(2.6) 



where 



i) (x, s) = <9 s lny9(x, s) . 



The problem (2.5), (2.6) is nonlinear. In addition both functions q and T are unknown here. 
Now the main question is How to approximate well both functions q and T using (2.5), (2.6)? 
Indeed, if we approximate them well (in a certain sense, specified below), then the target 
coefficient a(x) would be reconstructed easily via backwards calculations, see section 3. An 
equation similar with (2.5) was derived in the convexification method [6], [8]. However the 
major difference between our method and the convexification is in the numerical solution of 
the problem (2.5), (2.6). Indeed, it is solution of this problem which represents the major 
difficulty here. 

3 Layer Stripping With Respect to the Source Position 

An analogue of the nonlinear equation of this section for a different CIP, in which the original 
PDE was either hyperbolic or parabolic was previously derived in [3]. 

3.1 Nonlinear equation 

We approximate the function q (x, s) as a piecewise constant function with respect to the 
pseudo frequency s. That is, we assume that there exists a partition 



S = Sat < Sn_i < ... < Si < S = S, Sj_i — Si = h 




Vq (x, r) dr = (s n _i - s) Vq n (x) + h Vqj (x) , s E (s n , s n _i] . 



(3-1) 



We approximate the boundary condition (2.6) as a piecewise constant function, 

q n (x) =t(x),xG dQ, 



(3.2) 



where 



n-1 




(3.3) 



Hence, for s e [s n , s n -i) equation (2.5) can be rewritten as 



n-1 



L n (q n ) :=Aq n - 2Vq n ■ [h ^ Vq 3 - VT - 2h (Vq n f = 0. 



(3.4) 



In sections 4 and 5 we address the question on how to solve equations (3.4) for functions q n 
with the boundary conditions (3.2). 

3.2 Reconstruction of the target coefficient 

Suppose that functions {<?i}™ =1 are approximated via solving problems (3.2), (3.4) and that 
the tail function is also approximated. Then we reconstruct the target coefficient a (x) by 
backwards calculations as follows. First, we reconstruct the function v (x, s n ) by (2.4) as 

n 

v(x,s n ) = -h^2qj + T(x). (3.5) 

3=1 

In principle we can reconstruct the target coefficient a from (2.1). However, it is unstable to 
take second derivatives. Hence, we first reconstruct the function u (x, s n ) as 

u (x, s n ) = exp [v (x, s n )} . (3.6) 

Next, we use equation (1.1) in the weak form as 

— J VuVr) k dx = J aur] k dx, (3.7) 
n n 

where the test function r) k (x) , k — 1, K is a quadratic finite element of a computational 
mesh with r\ k (x) \qq— 0. The number K is finite and depends on the mesh we choose. 
Equalities (3.7) lead to a linear algebraic system which we solve. Then we obtain the function 

(au) (x) , 

K 

(au) (x) ^ 

k=i 

Hence, 



( x ) ~ Yl ak ^ ( x ) • 



a ( x ) ~ ~ / 1 - x akT lk ( x ) • ( 3 - 8 ) 



4 The Tail Function 



We consider in this section two procedures for obtaining sequential approximations for the 
tail function. First we find a first guess for the tail function using the asymptotic behavior 
of the solution of the problem (1.1), (1.2) as |x | — + oo, as well as boundary measurements. 
Second, we describe an iterative procedure with respect to tails. We call the combination 
of these two procedures "accelerators", because they help us to accelerate convergence of 
our method. We stress that we cannot prove convergence of the second procedure. However, 
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we have observed it in our numerical experiments. In our numerical experiments we have 
worked only with a rectangular domain. Hence, we assume in this section that 

(x, z) G Vl := {(x, z) : X\ < x < x 2 , z < z < B} . 

However, we do not yet know how to address the issue of tails in the case of an arbitrary 
convex domain Q. 



4.1 The first guess for tails 

First, we construct an approximation called "asymptotic tail". This is our first accelerator. 
We consider the fundamental solution of the problem (1.1), (1.2) for the case a(x,z) = k 2 . 
This solution is 

M = —K (k\(x -s,z- B)\), 

where K (z) a modified Bessel function. It is well known that the asymptotic behavior of 
this function is 



Ko(z) = i/^-^'(l + O(i), \z\ - oo. (4.1) 
V ^ \ z \ \ z \ 



Represent now solution of the problem (1.1), (1.2) as the solution of the following integral 
equation 

u(x,z,s) = —K (k\(x - s, z - z m )\) (4.2) 

-i- f K (k^/(x-0 2 + (z-v)^ [a V) - k 2 ] u V , s) d^d V . 
n 

Let S (x, z, s) = \(x,z) — (s,z m )\. The geometric meaning of S is illustrated in Figure 1. 
Introducing the function 

U(x, z, s) = 2V2nSe kS ■ u(x, z, s) 
and taking into account (4.1) and (4.2), we obtain that 
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l+g(x,z)+0[ ^ 



S — > oo. 



The function g(x, z) is unknown and is independent of S as S — > oo. Hence, we obtain for 
the function v — Inu 

v(x, z, s) = -kS+ l - ln(^) + g(x, z) + O (j^j , S -> oo, (4.3) 

where the unknown function g(x, z) is derived from g(x, z). 

We approximate the function g(x, z) by two different methods and the final answer is the 
average of two. The number of light sources N = 3 is taken in all our numerical experiments 
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when we approximate this function. We start at z = zq where the boundary values are 
known. We decompose the boundary values of v into 



1 7T 

v(x, zo, Sj) = -kSj + - In(^) + gj(x, z ) (4.4) 
for j=l,2,3. Then we average to obtain 

1 3 

g(x,z ,s) = - ^2gj(x,z ), (4.5) 

m=0 

Note that in (4.5) one should actually put "~" sign instead of "=". 

However, the above procedure (4.4)-(4.5) gives us the value of the tail functions v (x, z, s) 
only at z := z , i.e., v (x, z ,s) . Equation (4.3) provides an approximation for all (x, z) G Vl 
if we simply set g(x, z, s) = g(x, z , s). In our numerical experiments we found that this is 
insufficient. Hence, we use the measurement data from a different angle, which enhances our 
numerical results. We obtain a similar tail function using the measurement data at the lower 
edge of Q, i.e., at x — x\ and got a second tail function using the idea similar with the above. 
Thus, we have approximated v (xi,z,s) . Finally we set for the first guess for producing a 
tail function ^ 

Tx,o{x,z,s) := -[v{x,z Q ,s)+v{x 1 ,z,s)\. (4.6) 

4.2 The second accelerator: iterations with respect to tails 

The second accelerator involves another iterative process that enhances the reconstructed 
inclusion. Recall that k 2 := a is the constant background outside of our domain Q! . We 
now show how to find an approximation T\ (x, z, s) for the tail function. Let w^o = e v ° 
where Vq = T 10 (x, z, s) is the function introduced above. We reconstruct the approximation 
a± t i(x, z) for the unknown coefficient a(x, z) using the tail function (4.6) through the inversion 
formula in equation (3.7) for all quadratic finite element i] k : 



- J Vu lfi Vr] k dx = J a 1:1 u 1:0 r] k dx. 



Next, we apply (3.8). Then on the second step we solve the following boundary value problem 

Aii^i — ai t i (x, z) ui t i = 0, (x, z) G Q, 
ui,i \an= V (x, s) . 

The reason for doing so is that we need to satisfy the boundary condition obtained from 
measurements. 

We now describe a heuristic idea which motivates our iterative scheme. Let the function 
u be the solution of the following boundary value problem 

Au — a (x, z) u = —S(x, z), (x, z) G Vt, 
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u |an= s) 

with the unknown coefficient a (x, z) and the function -u satisfies 

Mo |an= 

with the background function a = k 2 . Denote p = u — u . Then 

Ap — a (x, z) p = (a (x, z) — a ) u , 
P \an= 0. 

Motivated by this idea, we introduce an iterative scheme and repeat the procedure until 
it converges. Suppose that after m — 1 iterations we have constructed the function Wi jm _i 
and have found the approximation ai tTn (x,z) > 0,m > 1 for the unknown coefficient a(x,z) 
using equation (3.7)-(3.8). Then on the iteration number m, we solve the following boundary 
value problem: 

Api, m - ai, m (x, z) p 1>m = (ai, m (x, z) - ai, m _i (x, z)) ui, m _i, 

Pi,m |an= 0. 

Next, we set 

To accelerate convergence, we modify the iterative scheme slightly to solve the following 
boundary value problems: 

Api, m - a hm (x, z) p hm = X m (a hm (x, z) - ai, m _i (x, z)) «i, m -i 

where 

exp{7r 2 e~ (m ~ 1) (ai, m (x, z) - a^m-i (x, z)) 2 } 

and 7 = 1.05. This choice of X m is made in numerical experiments. The choice of X m makes 
the sequence converge after about 50 iterations, instead of more than 300 in cases where 
X m = 1. 

Once we have u^ m = Wi >m _i +Pi, m , we construct ai ;7n+ i by equation (3.7) in the form of 



- J Vui, m V?7^x = J ai )m+ iiii >m 77,dx, 
n n 

for all quadratic finite elements = 1, , , , .,K and use (3.8) then. We iterate until the 
process converges, i.e., 

Il a l,mi — ft l,mi-l 1 1 X, 2 (f2) ^ 

\\ai, mi -i\\c(n) 



for a small e > of our choice, see (7.1) for a detail. We set for the first approximation for 
the unknown coefficient 

cii (x,z) := a 1>mi (x,z). 

Then we set for the tail 

Ti (x, z, s) = In u liTni (x, z) , (4.7) 

assuming that «i, mi > 0. Then we proceed with calculating the functions q n as in section 5. 

Remarks 4.1. 1. Unfortunately we cannot yet prove that functions a l m > 0. Therefore, 
we cannot prove analytically neither the existence of solutions of the above Dirichlet bound- 
ary value problems for functions pi >m nor the positivity of functions u 1}7n . Neither we cannot 
analytically prove that functions U\ jm converge, nor that our tail Ti is close to the correct 
tail T. Nevertheless, we observe all these "nice" properties in our computations. Figure 2 
displays the comparison of graphs of tails side by side, they have little visible difference. 2. 
Unlike [3], we do not change tails in all subsequent steps when calculating functions q n . In 
other words, the tail function is kept the same T := Ti (x, z,s) in all follow up steps of our 
algorithm. 



5 The Algorithm for Approximating Functions q n 

Step 1. We need to find an approximation for the function qi . To do this, we solve equation 
(3.4) with the boundary condition (3.2) at n — 1 iteratively for q 1 . That is, we should solve 

Aq 1 + 2VgiVTi = 2h (Vg x ) 2 (5.1) 

gi (x)=^(x),xGffi, (5.2) 
We solve the problem (5.1), (5.2) iteratively as 

Aq 1>k + 2Vgi )fc VTi - 2/iVgi )fc Vgi lfc _i = 0, q hk (x) = ^ (x) , x e dQ (5.3) 

where gi j0 = 0. 

We proceed with calculating the function (?i, m +i as in (5.3). We iterate in (5.3) until the 
process converges, i.e., 

||?i,fci - <7i,fci-iL 2 (Q) < e 

for a small e > of our choice, and e is the same as in (4.13). We set qi := qi tkl - The next 
reconstruction a 2 (x, z) is obtained using equations (3.5)-(3.8), where T := Ti. 

Step n. We now find an approximation for the function q n assuming that functions 
qi, ...,q n -i are found. We solve iteratively equation (3.4) with the boundary condition (3.2) 
as follows 

n-l 

Aq n , k - 2hJ^^(lj ■ Vg„, fc + 2Vg n , fe VT 1 - 2Wg„, fc Vg„, fc _i = 0,k = l, ...,m n , (5.4) 
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g n , fc (x) =t(x),xe dQ, (5.5) 
where g nj0 := g n _i.We iterate until the process converges, i.e., until 

||9n,fc n — Qn,k n -1 lli 2 (n) — £ 

for the above small e > 0. We set q n := q n ,k n - Then a n+ i(x, z) is reconstructed using 
equations (3.5)-(3.8), T := 7\. 

We find functions oi, a at for n = 1, N, where N is the number of sub intervals of the 
interval [so,s] . Finally, the resulting function a (x, z) is 

1 N 

a(x,z) = —^2ai(x,z) . (5.6) 
i=i 

We stress that we did not prove convergence of tails T n nor q n rigorously. Neither we 
cannot prove that functions p ntm in section 4 are positive, because we cannot prove that 
a>n,m — a>n,m-i < (in order to apply the maximum principle). However, we have observed 
both the positivity of functions p n , m and convergence of tails T n and q n in our computations. 



6 Convergence 

Below we follow the concept of Tikhonov for ill-posed problems [10], which is one of backbones 
of this theory. By this concept one should assume first that there exists an "ideal" exact 
solution of the problem with the exact data. Next, one should assume the presence of an error 
in the data of the level (, where ( > is a small parameter. Suppose that an approximate 
solution is constructed for each sufficiently small (. This solution is called a "regularized 
solution", if the (— dependent family of these solutions tends to that exact solution as q 
tends to zero. Hence, one should prove this convergence (Theorem 6.1). 

In this section we use the Schauder's theorem to estimate functions q n .k-, see §1 of Chapter 
3 of [7] for this theorem. Since the Schauder's theorem requires C 2+a smoothness of the 
boundary dQ, we assume in this section that Q C 1R 2 is a convex bounded domain with 
dQ G C 2+a . This is, of course in a disagreement with the above case of Q being a rectangle. 
However, we use the rectangle only because of the problem with tails, in which we cannot 
rigorously prove that they are small and do not yet know how to approximate them well 
heruistically for the case of a more general domain, so that they would be close to correct tails. 
However, an analogue of our convergence result (Theorem 6.1) can be proven for the case 
when Q is rectangle and an FEM (i.e., discrete) version of equation (3.4) is considered with a 
fixed number R of finite elements. To do this, one needs to consider the weak formulation of 
(3.4) and to use the Lax-Milgram theorem instead of the Schauder's theorem. Although the 
Lax-Milgram theorem would provide only estimates of H 1 norms of functions q n rather than 
more desirable C 2 norms, but using the equivalency of norms in finite dimensional spaces, 
we can still get estimates of C 2 norms and these estimates would naturally depend on R. 
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6.1 Exact solution 

Following the Tikhonov concept, we need to introduce the definitions of the exact solution 
first. We assume that there exists an exact coefficient function a* (x) G C a (fi) , a = const, G 
(0, 1) , which is a solution of our Inverse Problem. Let the function 

u* (x,s) gC 2+q? (|x-x | > e),Ve > 0,Vx = (s, B) > 0, Vs G [s,s] 

be the solution of the problem (1.1), (1.2) with a (x) := a* (x). Let 

v* (x, s) = In m* (x, s) , g* (x, s) = — ^ ^ , T* (x, s) = (x, s) . 

By (2.1) 

a* (x) = Av* + (Vv*) 2 . (6.1) 
Also, the function q* satisfies the following analogue of equation (2.5) 

s 

Aq* - 2V<f • J V<f (x, r) dr + 2V<fVT* = (6.2) 

s 

with the boundary condition (see (2.6)) 

q* (x, s) = ip* (x, s) , (x, s) G <9f2 x [s, s] , (6.3) 

where ip* (x, s) = <9 S In y?* (x, s) , where </?* (x, s) = u* (x, s) for (x, s) G dQ x [s, s] . 

Definition. We call the function q* (x, s) the exact solution of the problem (2.5), (2.6) 
with the exact boundary condition ip* (x, s). Naturally, the function a* (x) from (6.1) is called 
the exact solution of our Inverse Problem. 

Therefore, 

q* (x, s) G C 2+a (H) x C 1 [s, s] . (6.4) 

We now approximate the function q* (x, s) via a piecewise constant function with respect to 
s G [s, s] . Let 

Sn— 1 ^n — 1 

9n ( x ) = \ J Q* ( x ' s ) ds > C ( x ) = \ J ^* (^ s ) ds 

Then by (6.4) 

q* (x, s) = q* (x) + Q n (x, s) , ip* (x, s) = ijj* n (x) + ^ n (x, s) , (6.5) 
s G [s n , s n _i] , where functions Q„, ^„ are such that for s G [s n , s„_i] 

ll<5n(x,s)|| c2+a ^ < C*/i, ||* n (x,s)|| c2+a ^ <C*/i,VsG [s„,s„_i] ,ra = 1, ...,iV, (6.6) 



12 



where the constant C* > depends only on C 2+a (Q) x C 1 [s, s] and C 2+a {dVt) x C 1 [s, s] 
norms of functions q* and ip* respectively. Hence 

<£(x)=<(x),xeafi (6.7) 

and the following analogue of equation (3.4) holds 

L n (q n ) := Aq* n - 2Vq* n ■ £ Vg* - VT*j - 2h {Vq* n f = F„* (x, h) . (6.8) 



where the function F n (x, h) G C a (fi) and 

\\K(^h)\\ Ca{ ^<C*h. (6.9) 

We also assume that the data ip (x, s) in (1.3) are given with an error. This naturally 
produces an error in the function tp (x, s) in (2.6). An additional error is introduced due to 
taking the average value of if) (x, s) over the interval (s n ,s„_i). Hence, it is reasonable to 
assume that 

|C (x) - i, n (x)|| c2+Q(m) <Ci(a + h), (6.10) 

where a > is a small parameter characterizing the level of the error in the data ip (x, s) 
and the constant C\ > is independent on numbers a, h and n. 

Remark 6.1. It should be noted that usually the data ip{x,s) in (1.3) are given with 
a random noise. Although the differentiation of the noisy data is an ill-posed problem, but 
there exist effective numerical regularization methods of its solution. We are not addressing 
the corresponding theory here referring the reader to e.g., [4], and also see section 7 for our 
way of handling it. 

6.2 Convergence theorem 

First, we reformulate the Schauder's theorem in a way, which is convenient for our case. 
Introduce the positive constant M* as 



AT 



max 

Kn<7V 



(ll?nllci+«(n)) + 2 W T *Wc^(n) + H V <7illc«(n) + 1 



, C*, Ci > , 



where C* and C\ are constants from (6.9) and (6.10) respectively. Consider the Dirichlet 
boundary value problem 

3 

Au + bj{x)u Xj — d{x)u = f (x) , x £ f2, 

3=1 

u\ dn =g(x),geC 2+a (dQ), 
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where functions 



b v d,feC a (Q),d(x)>0; max(||y ca(n) Jd|| ca(n) ) < 4M*. 



By the Schauder theorem there exists unique solution u G C 2+a (fl) of this problem and 
with a constant K = K (M*, Vl) > the following estimate holds 



\9\\c 2 + a (dn) + 11/ lie 



In Theorem 6.1 we use a function T appr (x, z,s) instead of the above constructed function 
T\ (x, z,s) only because the latter was constructed for a rectangle, while Theorem 6.1 works 
with a convex bounded domain, also see the beginning of this section. 

Theorem 6.1. Let Q C 1R 2 be a convex bounded domain with the boundary dfl G C 3 . 
Suppose that an approximation T appr (x, z, s) for the tail is constructed in such a way that 

\\T appr — T*\\ C 2 +a ^ < £, (6-11) 

where £ G (0, 1) is a sufficiently small number and that this function T appr (x, z,s) is used in 
(5.3), (5.4) instead of the function Ti (x, z,s) . Denote r] = h + a + ^ + e. Suppose that the 
number f3 := s — s = Nh is such that 

B< . (6.12) 

Then there exists a sufficiently small number i] = r) (K (M*,Q) ,M*,c,s,s) G (0,1) and 
a sufficiently large small number h = h (K (M* , Q) , M* , c, s,s) G (0,1) such that for all 
i] G (0, r] ) and for every integer n G [1, N] the following estimates hold 

Ik - q:\\ c ^(n) < 2KM * ( h + > ( 6 - 13 ) 
\\q n \\ c i +a{ n) <2M*. (6.14) 

Remark 6.2. As it was stated above, unlike the time dependent case of [3], we cannot 
prove the estimate (6.11) for tails. However, we observe convergence of tails in computations 
if taking T := T\ as in (4.7). 

6.3 Proof of Theorem 6.1 

In the course of this proof we assume that 77 G (0, rj ) , h G (0, r] ) . Denote 

5n,fe(x) = g„,fc(x) - g*(x),f (x) = T appr (x) - T*(x),^ n = $ n - ip* n , (6.16) 

n-l 

v n , k (x, s n ) = -hq n , k -h^2qj + T appr (x) , u n<k (x, s n ) = exp [u n>fc (x, s n )] , (6.17) 
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v n , k (x, s n ) = v n , k (x, s„) - < (x, s n ) , a n , fc (x) = a nk (x) - a* (x) . (6.18) 

The proof basically consists in estimating these differences. 

First, we estimate q ltl . Set in (6.8) n = 1 and subtract it from (5.3) at k — 1, recalling 
that q lfi = 0. We obtain 

Agi,! - 2Vq 1A VT appr = 2V gi *VT - 2h (Vq*) 2 - F*, 

|an= ^i- 

By Schauder theorem, and (6.9)-(6.11) we obtain 

||5i,i|lc^(n) <^M*(/i + 77). (6.19) 

Hence, 

||?i,i + Silicon) = IIMIc^n) < M* + XM* (/i + t?) < 2M*. (6.20) 
Now we estimate q-y^. Set in (6.8) n — 1 and subtract it from (5.3). Then 

Agi >fc + 2 (VT appr - hVq 1>k -i) Vgi, fc = (6.21) 
-2V 9l *(vf -2Wgi ifc _i) -F*, 

and also 

|an= (6-22) 
Set in (6.21) k — 2. Then (6.19) and (6.20) imply that 

||2 (VT appr - hVq ltl )\\ Ca( fl < 2 (M* + 2/iM*) < 3M*, (6.23) 

2Vg x * fvf - 2hVqi A ) + < 2AT [f + 2#M*/i (/i + 77) + rj/2] . 

By (6.12) 2KM*h < 1/2. Hence, 

2Vq{ fvf - 2hVqi 1) + ,_ x < 2M* (/i + 2?j) . 

Hence, by (6.21)-(6.23) and Schauder's theorem 

fellc^jj) <2KM* (h + 3r]). 

and similarly with (6.19) 

\\qiM\c^(n)< 2M *- 

Assume that 

\\qi,k-i\\ c2+a (n) < 2KM* (h + 3 V ) , ||gi >ifc -i|| c2+a(n) < 2M*. (6.24) 
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We now estimate the function qi tk -i- Similarly with the above 

||2 (VT appr - hVq 1 , k ^)\\ Ca(Ti) < 2 (M* + 2hM*) < QM*. 

Next, using (6.24), we obtain 

2Vql fvf - 2hVq hk - 1 ) + F* < 2M* [f + 2/TAT/i (/i + 77) + rj/2] . 

V / C a [€i\ 

By (6.12) 2fTM*/i {h + rj) < 1/2 (/i + 77) . Hence, 

£ + 2KM*h (h + 77) + 77/2 < /i + 277. 

Hence, 



(6.25) 



< 2M* (/i + 277) . 



(6.26) 



2Vq{ (VT - 2hVq 1 , k - 1 ) + F* 

Hence, Schauder's theorem, (6.21), (6.25) and (6.26) lead to 

\\qi,k\\c^(n) <2KM*(h + 3r l ),\\q 1 , k - 1 \\ ca+a ^ < 2M*, k = 1, 2, ... (6.27) 

We now estimate the function q n>k , assuming that (6.27) holds for functions cfi,qi with 
j < n, as well as for functions q n ,m-> qn.m with m < k — 1. In other words, we assume that 



hWcz+^n) < 2KM * (h + 3r/) , H\\ c2+a( n) < 2M*. 



and 



1 1 q n ,m 1 1 C 2+«(n) ^ 2fTM*( /1 + 377), I |g n>m || C 2+ a (n) < 2M*,m < A; - 1. 
Subtracting (6.8) from (5.4) and using (5.5), (6.3) and (6.16), we obtain 



ra-l 



Mn.k - 2 ( h^Vqj - {VT appr - 2hVq n ,k-i) ) Vq n . k 
3=1 



(6.28) 
(6.29) 

(6.30) 



n— 1 



2 ^ v% v g ; + 2V g ;VT + 2w g ;v? n ,_ 1 , 



Estimate first the coefficient at Vq n . k in (6.30). Using (6.28) and (6.29), we obtain 

t-i 



(6.31) 



hJ2^q 3 



3=1 



< AM*Nh, 



2 ||VT appr - 2W ?n , fe _ 1 || Ca(n) < 2 (M * + 2hM*) < 3M*. 
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Since by (6.12) AM*Nh < AT, then the estimate for that coefficient is 



n— 1 



+ 2 ||VT appr - 2/iVg n , fc _i|| Ca ( ? j) < 4AT. 



(6.32) 



Hence, we can apply Schauder's theorem with the constant K. Now we estimate the right 
hand side of equation (6.30). Using (6.28) and (6.29), we obtain 

n-l 



< AK (M*) Nh (h + 3r)) , 



||2W<£Vg n . fc _i|| Ca(n) < (M*) 2 h(h + 3rj) , 
2Vq* n VT / x < 2ATr/. 
Estimate the right hand sides of the last three inequalities. By (6.12) we have 



4K (AT) 2 Nh (h + 3?7) < ^- (h + 3??) , 

9 AT 
4K (M*) h(h + 3 V )<—(h + 377) . 



M* 



Hence, 



n— 1 
J'=l 



c«(n) 



+ I^W^V^II^ + phVq'Vqn.k-^p) (6-33) 



M* ( h 3 

< — (/i + 3*7) + 2M* V = M* (- + -r)) . 



By (6.30)-(6.33) and Schauder theorem we obtain 

\\qr,k\\ c ^(n) < KM * Q + ^) + K V< KM* Q + ^ < KM* (h + 3 V ) . (6.34) 
Hence, 



\\q n .k\\c2+ a (Ti) = hn.k + q*Jc*+<*(n) < KM * ( h + 3r ?) + AT < 2AP. 
Estimates (6.34) and (6.35) complete the proof of this theorem. □ 



(6.35) 



7 Numerical Studies 

We have performed numerical experiments on several cases of reconstructions using the 
method discussed above. We have chosen the range of geometrical parameters of the rectangle 
n, which is typical for optical imaging of small animals and have chosen the range of optical 
parameters typical for biological tissues [1] ,[11], [13] . 
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7.1 Some Details of Numerical Studies 



For the forward problem, we calculate the solution of the diffusion equation 

DAu — fi a (x, z)u = —5 (x — s, z — z m ) (7.1) 
with the conventional condition at the infinity 

lim u(x, z, s) — 0, (7.2) 

\(X,Z)\-HX> 

where D — 1/ (3//J = const. > is the diffusion coefficient, where optical coefficients fi' s 
and /i a (x,z) were discussed in Introduction (Section 1). In our computations the function 
fi a (x, z) is unknown and the constant ji' s is given. Thus, in our case a(x, z) = 3fi' s ■ fi a (x, z) 
(compare with (1.4)). Consider the rectangle Q, 

Q = {(x, z) : 5cm < x < 15cm, 5cm < z < 10cm} . 

We assume that 

a(x, z) = k 2 = const. > in M 2 \a (7.3) 

We assume that in (7.1) the source position (s, z m ) is running along the right side of Q, i.e., 
z m = L = 10cm. Also, consider a bigger rectangle 

Q = {(x, z) : 0cm < x < 20cm, 0cm < z < 15cm}. 

The reason why we consider the rectangle Qq along with the rectangle Q is that it is natural to 
approximate the solution of the problem (7.1), (7.2) in the infinite domain by the solution of 
equation (7.1) in Qq with Robin boundary conditions at dQo. We have established numerically 
that for the range of parameters we use the solution of the problem (7.1), (7.2) is close in fl to 
the solution of equation (7.1) in the bigger rectangle fl with the Robin boundary conditions 
at its sides. Figure 3 illustrates rectangles Q and Q. 

The light sources are located in several positions (xi,z) = (si, 10) along the right-hand 
side of the smaller rectangle Q, and receivers, which mimic the so-called CCD camera are 
located at the left-hand side of Q. CCD stands for a "charge-coupled device". A CCD 
camera is an image sensor, consisting of an integrated circuit containing an array of linked, 
or coupled, light-sensitive capacitors. A typical CCD camera can take up to 512 x 512 data 
points simultaneously, which will provide an adequate amount of data for our reconstruction. 
In all three examples, we have used an ideal light source modeled by the function —5 (x — 
Si,z — 10) in the 2D case of (1.1). In numerical simulation 5(x — si,z — 10) = cn(si, 10), 
where n is the finite element at the location, and c is the scaling constant to ensure that the 
area equals one. 

We use three (3) sources to construct an approximation of the tail functions which was 
described above. Next, we use all five (5) sources for the above layer stripping procedure 
both in the s-derivative and the s-integral. 
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We have generated the data for the forward problem for total of five (5) different locations 
of the light source, Sj = 1, 5, where si = 0, Sj = Sj_i + 0.625cm, % = 2, 5. Hence, N = 4 
and we have used four (4) functions q n . An increase of the number N did not result in 
significant improvements of results. Note that a similar observation took place in numerical 
experiments of [3] for the case of an increase of the number of functions q n after a certain 
"limit" . In our reconstruction method, we use the solution of the forward problem to generate 
the data for the inverse, add noise to the measurement data, and reconstruct the absorption 
coefficient /i a (x,z) in Q. The domain Q will be our basic computational domain for our 
inverse calculations. 

In our examples, the coefficients in equation (7.1) are D = 0.02cm uniformly and 
fi a = 0.1cm -1 at all grids except off the inclusions, and in inclusions ji a ranges from 0.1 
to 0.3 cvrT 1 . The maximum inclusion/background contrast is 3:1 in our computations. Our 
algorithm calculates the forward problem with Robin boundary conditions at 8Qq, given the 
distribution of the absorption coefficient. A total of 130x93 rectangular finite elements is 
used for forward calculations for the domain Q - 

For the simulated boundary measurements we take the solution of the forward problems 
along the left and lower boundaries of Q to construct first guess for tails (section 4.1) and 
then we take all 4 sides for the rest of the problem. The number of measuring points is 65 
along the left edge of Q and 31 along the lower edge of Q. The number of measuring points 
at the low left corner is shared by both sides and therefore the total number of independent 
measuring points is 95. 

For each detector position, we introduce the random noise as the random process with 
respect to the detector locations, y3(x, Sk) = y?(x, Sk) [1 + x ( x )] > where x ( x ) is the random 
variable, which we introduce as x — 0.02W, where W is a white noise with equal distribution 
at [-1,1]. Hence, this is 2% of the multiplicative random noise. To obtain the realistic first 
s-derivatives, we started with simulated light distribution u(x, z, s) added with similar noise 
to simulate the situation used in applications and let v = lnw. Then we take first derivatives 
with respect to s as shown in the paragraph below. 

A regularization method was introduced to pre-process the noise in the measurement 
data. We use a polynomial approximation with respect to the detector location x. In 
our setting, the measurements are collected at 65 locations along the left boundary and 31 
points along the lower boundary. We use an eight order polynomial to approximate functions 
</3(x, Sk) with respect to x for each s*.. The polynomial is optimal in the least square sense 
[12], [14], and its sub-routing is commonly available, see for example 

http://perso.orange.fr/jean-pierre.moreau/fJstsqr.html. We demonstrate the essence of 
the approximation in Figure 4. Thus, we have obtained approximate polynomial functions 
Tp(x,Sk). We use functions Tp(x,Sk) instead of ip(x,Sk). The first s-derivatives are processed 
afterwards by the formula 

f( Sl ) « /(g2) ~ /(gl) . 

*2 — s l 

and similar ones for other source locations. 
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7.2 Numerical Experiments 

In the following numerical examples, we illustrate the results in a few different shapes and 
locations of the two inclusions. Our method has shown its success in dealing with those 
cases. In all cases, iterations with respect to tail (as described in section 4.2) were only done 
for the first light sources. Additional light sources did not bring any significant changes to 
reconstruction, therefore iterations with respect to tails are not shown in this paper. The 
total number of elements K in (3.8) is 450 in our calculation. 

The Convergence Criterion for functions a m in the procedure of finding the second accel- 
erator Ti for tails is 



\a lmi (x,z) 



. E 



(^l,mi %j) ^l,mi-l %j) ) I 



■ ,Jmax 



/ 3^max|ai )mi _i(a; i ,2!j 



(7.1) 



where Ni = i ma xjmax is the total number of finite elements. In all our examples, e = 10' 
The number of iterations required for convergence is listed below: 

The number of iterations required for convergence 



Example 1 


Example 2 


Example 3 


52 


50 


48 



Example 1. Inclusions are two circles with the radius 1 cm, and their centers are placed 2 
cm off the left edge. The coefficient is fJ, a (x, z) = 0.3 inside inclusion and /i a (x, z) = 0.1 = k 2 
outside of inclusions. We have also added 2% of random noise to the boundary measurements, 
see subsection 7.1. 

Figure 5a displays the original distribution and its 1-d cross section. Figure 5b shows 
reconstruction from the noisy data and its 1-d cross section. 
The relative errors of the reconstruction are as follows: 
Table 1. The relative errors of reconstructions in Example 1 



RMSE 


AME 


ME 


0.312366805537619 


0.115817263155519 


-0.058254402556204 



Note that for the data set (xi, X2, • • • , xjvj and its approximation (xi,X2, • • • , xni), the 
values of the function /i a (x, z) taken at each of the grid points, the Relative Root Mean Square 
Error (RMSE), Relative Absolute Mean Error (MAE) and Relative Mean Error (ME) are 
calculated by 



RMSE = 



Y.{x k -x k ) 2 J2\x k 
MAE = k=l 



x k \ 



r Nl max \xk\ 



N 1 max \Xk 



-, ME = 



k=l 

Ni max \xk\ 



In our case x k are correct values of the coefficient fi a (x, z) at the grid points of the 
sub-rectangle 

Q! = {(x, z) : 5cm < x < 15cm, 5cm < z < 8cm} C Q. 
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We illustrate in Figure 5c the difference of two consecutive reconstruction 



\a m (x,z) - a ro _i(a;,2)|| = 



ViVT max | a m _i(xi,^-) | 

as a function of the number of iteration to. Figure 5d depicts the relative error in comparison 
with actual inclusion expressed by 



Y, \(a m (xi, Zj) - a(xi, Zj))\ 2 

RMSE = -? — — 

V^i max |a(a;i, Zj) | 

as a function of the number of iterations to. 

Example 2. Inclusions are two circles of the radius 1cm, and their centers are placed 2 
cm off the left edge. The coefficient fi a (x, z) is defined as 

i \ — \ max [0.3 cosd(x, z), 0.1] , inside of each circle 1 , ■. 

^ x i z > = \ ^ o.l otherwise, J [ ' 

where d(x, z) is the minimum distance to center of each of these two circles. 

Figure 6a displays the original function in two inclusions and its 1-d cross section. Figure 
6b shows the reconstruction result with 2% noise and its 1-d cross section. 

The relative errors of reconstruction are as follows: 

Table 2. The relative errors of reconstruction in Example 2 



RMSE 


AME 


ME 


0.261946376827213 


0.087514613280764 


-0.024755081675791 



We illustrate in Figure 6c the difference of two consecutive reconstruction \ \a m (x,z) — 
a m -i(x, z)\ \ as a function of the number of iteration m. Figure 6d depicts the relative error 
in comparison with actual inclusion expressed by RMSE as a function of the number of 
iterations m. 

Example 3. Inclusions are two circles of the radius 1 cm and 0.6 cm, whose centers are 
placed 2 cm off the left edge. The coefficient [i a (x, z) is defined as 

( \ — \ max [0-3(cos d(x, z)(l + 0.1r7(a:, z)), 0.1] , inside of each circle 1 , . 
^ z ) = \ 0.1 otherwise. J (7 ' 5j 

Similarly with (7.4) d(x, z) is the distance to center of the circle. In (7.5) rj is a realization 
of a white noise valued between [-1, 1]. The random pattern is introduced to test the ability 
of our method to handle complex shapes. See Figures 7a,b for results. 
Table 3. The relative errors of reconstruction in Example 3 



RMSE 


AME 


ME 


0.336142461513025 


0.084080392591033 


-0.026699370402152 



We illustrate in Figure 7c the difference of two consecutive reconstruction \ \a m (x,z) — 
a m -i(x, z)\ \ as a function of the number of iteration m. Figure 7d depicts the relative error 
in comparison with actual inclusion expressed by RMSE as a function of the number of 
iterations to. 
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Figure 1. The figure illustrates construction of the real distance S in the asymptotic 
expansion of the tail function. The distance S measures the distance to any interior point 
individually, rather than the distance to the edge denoted by s. 




■ :%T.S 

■L.!.-1cJ 
'22^1 

.67825 
<IC;> /2 



Figure 2. These two figures illustrate the comparison of actual tail function T (x, z, s) 
from forward problem (left panel) and the calculated tail function T\ (x, z, s) derived from 
iterative procedure (right panel). 
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Figure 3. The figure illustrates relative positions of rectangles Vt$ and Vt. The light source 
location is at the right edge Q. 




Figure 4- We use an eighth order polynomial to approximate functions <p(x,Sk) and 
ip(x,Sk) with respect to x for each Sk- The polynomial is optimal in the least square sense, 
and its sub-routing is commonly available [ 11,12 }. We demonstrate the essence of the 
approximation in Figure 4- 
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Iteration of Improving Tail 



Iteration of Improving Tail 



Figures 5a-5d. Figure 5a displays the original distribution. Figure 5b shows reconstruc- 
tion from the noisy data. The 1-d cross sections next to the 2-d figures show the profiles of 
the original inclusion and its reconstruction at z = 7 

We illustrate on Figure 5c the difference of two consecutive reconstruction 



\{a m (xi, zj) - a m -i(xi, Zj))\ 2 

i = 1 j ■ ■ ■ )1>Tnaxij = lj ■ ■ ■ ijmax 



r N\ max \ a m -i(xi, z 
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as a function of the number of iteration m. The function is used for determining stopping 
criterion. Figure 5d depicts the relative error in comparison with actual inclusion expressed 
by 

J E _ \(am(xi,Zj) - a(x h Zj))\ 2 
RMSE = -I 



V^Vi max \a(xi, Zj) 
as a function of the number of iterations m. 
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Iteration of Improving Tall 



Figures 6a-6d. Figure 6a displays the original function in two inclusions. Figure 6b 
shows the reconstruction result with 2% noise. We illustrate on Figure 6c the difference of 
two consecutive reconstruction 



\(a m (xi,Zj) - a m -i(xi, Zj))\< 

i— 1 5 ... )imax ij ' — 1 1 ■ ■ ■ ijmax 



r Nxmax.\a m -i(x i , Zj) 
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as a function of the number of iteration m. The function is used for determining stopping 
criterion. Figure 6d depicts the relative error in comparison with actual inclusion expressed 
by 

J E _ \(am(xi,Zj) - a(x h Zj))\ 2 
RMSE = -I 



V^Vi max \a(xi, Zj) 
as a function of the number of iterations m. 
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Figures la-Id. Figure 6a displays the original distribution. Figure 7b displays the recon- 
struction result with 2% noise in the data. The 1-d cross sections show the profiles of the 
original inclusion and its reconstruction at z = 7 We illustrate on Figure 7c the difference 
of two consecutive reconstruction 



J2 \{a m (xi, zj) - a m -i(xi, Zj))\ 2 

i = 1 j ■ ■ ■ )imax ij — 1 5 ■ ■ ■ ijmax 



r N\ max \ a m -i(xi, z 
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as a function of the number of iteration m. The function is used for determining stopping 
criterion. Figure Id depicts the relative error in comparison with actual inclusion expressed 
by 

J E _ \(am(xi,Zj) - a(x h Zj))\ 2 
RMSE = -I 



V^Vi max \a(xi, Zj) 
as a function of the number of iterations m. 
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